4  Psykiska/psykosomatiska besvär

Author
Affiliation
Magnus Johansson
Published

March 13, 2023

4.1 Bakgrund

Item/frågor har etiketter F88-F99 i datafilen, och motsvaras av fråga 90-101 i PDF-filen med frågor.

Samtliga frågor har fem svarskategorier, vilka varierar mellan frågorna. Fem frågor har svarskategorier från “Aldrig” till “Flera gånger i veckan”. Sex frågor har från “Sällan” till “Väldigt ofta”, och en från “Inte alls” till “Väldigt mycket”.

Svarsdata har kodats så att högre poäng innebär mera besvär/högre risk.

Sektionen i enkäten inleds med meningen: “NÅGRA FRÅGOR OM HUR DU MÅR”.

4.2 Lista med enkätfrågorna

Code
RIlistitems(df.omit.na)
itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F89 Känner du dig ledsen och deppig utan att veta varför?
F90 Händer det att du känner dig rädd utan att veta varför?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94 Hur ofta tycker du att du inget duger till?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F96 Är du nöjd med ditt utseende?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?

4.2.1 Demografi

Vi har 9595 deltagare i samplet från 2014, och deras könsfördelning återges i tabellen nedan. Deltagare som saknar data på samtliga frågor är borttagna ur analysen.

4.3 Deskriptiva data

4.3.1 Demografi

Code
RIdemographics(dif.gender, "Kön")
RIdemographics(dif.arskurs, "Årskurs")
Kön n Percent
Flicka 5179 54
Pojke 4416 46
Årskurs n Percent
Åk 9 4109 42.8
Gy 2 5486 57.2

4.3.2 Item-data

Code
RItileplot(df.omit.na)

Code
RIbarstack(df.omit.na)

Code
RIbarplot(df.omit.na)

4.4 Rasch-analys 1 samtliga items

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F89 Känner du dig ledsen och deppig utan att veta varför?
F90 Händer det att du känner dig rädd utan att veta varför?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94 Hur ofta tycker du att du inget duger till?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F96 Är du nöjd med ditt utseende?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIitemfitPCM2(df.omit.na, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F88 1.06 1.043 0.749 0.775
F89 0.709 0.716 -4.108 -4.521
F90 0.851 0.896 -0.773 -1.232
F91 1.133 1.094 1.129 1.354
F92 0.788 0.794 -2.92 -3.126
F93 1.029 1.015 0.405 0.38
F94 0.839 0.819 -2.234 -2.329
F95 1.059 1.054 0.924 0.694
F96 0.995 0.969 0.27 -0.508
F97 0.968 0.959 -0.389 -0.415
F98 0.956 0.965 -0.557 -0.561
F99 1.072 0.995 0.507 -0.099
Code
RIpcmPCA(df.omit.na)
PCA of Rasch model residuals
Eigenvalues
2.02
1.31
1.22
1.11
1.08
Code
RIloadLoc(df.omit.na)

Code
RIresidcorr(df.omit.na, cutoff = 0.2)
F88 F89 F90 F91 F92 F93 F94 F95 F96 F97 F98 F99
F88
F89 -0.09
F90 -0.06 0.11
F91 -0.05 -0.13 -0.1
F92 -0.2 -0.01 -0.08 -0.12
F93 0.02 -0.08 -0.06 0 -0.1
F94 -0.18 0.06 -0.03 -0.21 0.12 -0.15
F95 -0.02 -0.17 -0.11 -0.02 -0.2 -0.06 -0.24
F96 -0.19 -0.12 -0.12 -0.22 0.29 -0.21 0.17 -0.24
F97 -0.11 -0.02 -0.05 -0.16 -0.07 -0.17 0.02 -0.13 -0.01
F98 0.02 -0.11 -0.04 -0.07 -0.21 -0.06 -0.17 0.18 -0.22 -0.14
F99 -0.19 0.02 -0.1 -0.13 -0.01 -0.17 -0.03 -0.13 0.08 -0.02 -0.17
Note:
Relative cut-off value (highlighted in red) is 0.122, which is 0.2 above the average correlation.
Code
RItargeting(df.omit.na)

Code
RIitemHierarchy(df.omit.na)

Code
plot(mirt.rasch, type="trace")

PCA på residualer från Rasch-modellen indikerar möjlig multidimensionalitet. Utifrån faktorladdningarna på första residualkontrasten ser vi två potentiella kluster i data. Det ena betecknas av de psykosomatiska frågorna, det andra av psykiska besvär. Vi tittar på klustren separat.

4.5 Rasch-analys 1 - psykosomatiska besvär

Code
items.psm <- c("F88","F91","F93","F95","F98")

df.psm <- df.omit.na %>% 
  select(any_of(items.psm))

RIlistItemsMargin(df.psm, fontsize = 13)
itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
Code
RIitemfitPCM2(df.psm, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F88 0.897 0.894 -1.126 -1.216
F91 0.908 0.919 -0.996 -0.98
F93 0.894 0.889 -1.127 -1.712
F95 0.783 0.802 -2.589 -2.966
F98 0.77 0.79 -3.037 -2.947
Code
RIpcmPCA(df.psm)
PCA of Rasch model residuals
Eigenvalues
1.50
1.28
1.17
1.04
0.01
Code
RIresidcorr(df.psm, cutoff = 0.2)
F88 F91 F93 F95 F98
F88
F91 -0.23
F93 -0.13 -0.17
F95 -0.22 -0.24 -0.26
F98 -0.14 -0.28 -0.23 0.02
Note:
Relative cut-off value (highlighted in red) is 0.012, which is 0.2 above the average correlation.
Code
RItargeting(df.psm)

Code
RIitemHierarchy(df.psm)

Code
plot(mirt.rasch, type="trace")

Code
for (i in 1:ncol(df.psm)) {
    barplot(table(df.psm[, i]), col = "#8dc8c7", main = names(df.psm[i]), 
      ylab = "Number of responses", xlab = (itemlabels %>% 
                                              filter(itemnr %in% names(df.psm)) 
                                            %>% .[i,2]))
  }

Bortsett från svarskategorierna fungerar items acceptabelt.

4.5.1 Omkodning av svarskategorier

Vi behöver åtgärda flera items genom att slå samman svarskategorier:

  • F91: 1+2 och 3+4
  • F93 och F98: 3+4
  • F95: 0+1
Code
df.psm$F91 <- recode(df.psm$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.psm$F93 <- recode(df.psm$F93, "4=3", as.factor = FALSE)
df.psm$F98 <- recode(df.psm$F98, "4=3", as.factor = FALSE)
df.psm$F95 <- recode(df.psm$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)
Code
plot(mirt.rasch, type="trace")

Code
RIitemHierarchy(df.psm)

4.6 Rasch-analys 2 - psykosomatiska besvär

Efter att ha åtgärdat svarskategorierna.

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
Code
RIitemfitPCM2(df.psm, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F88 0.886 0.891 -1.438 -1.2
F91 0.937 0.939 -0.624 -1.136
F93 0.904 0.901 -0.925 -1.528
F95 0.832 0.84 -2.101 -2.423
F98 0.798 0.816 -2.623 -2.55
Code
RIpcmPCA(df.psm)
PCA of Rasch model residuals
Eigenvalues
1.47
1.30
1.14
1.09
0.00
Code
RIloadLoc(df.psm)

Code
RIresidcorr(df.psm, cutoff = 0.2)
F88 F91 F93 F95 F98
F88
F91 -0.18
F93 -0.19 -0.08
F95 -0.27 -0.14 -0.28
F98 -0.23 -0.18 -0.26 -0.05
Note:
Relative cut-off value (highlighted in red) is 0.013, which is 0.2 above the average correlation.
Code
RItargeting(df.psm)

Code
RIitemHierarchy(df.psm)

4.7 Invarians-/DIF

4.7.1 Kön

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
Code
RIdifTable(df.psm, dif.gender)

Item 2 3 Mean location StDev MaxDiff
F88 -0.107 0.033 -0.037 0.099 0.140
F91 0.165 0.131 0.148 0.024 0.034
F93 -0.216 -0.140 -0.178 0.054 0.076
F95 0.103 -0.130 -0.014 0.164 0.232
F98 0.056 0.105 0.080 0.035 0.050
Code
RIdifFigure(df.psm, dif.gender)

Inga problem med DIF för något item gällande kön.

4.7.2 Årskurs

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
Code
RIdifTable(df.psm, dif.arskurs)

Item 2 3 Mean location StDev MaxDiff
F88 -0.088 0.005 -0.041 0.066 0.093
F91 0.079 0.197 0.138 0.084 0.118
F93 -0.140 -0.217 -0.179 0.055 0.077
F95 -0.017 -0.007 -0.012 0.007 0.010
F98 0.166 0.022 0.094 0.102 0.145
Code
RIdifFigure(df.psm, dif.arskurs)

4.7.3 Årtal

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
Code
final.psm.items <- names(df.psm)
# write.csv(final.psm.items, file = "2022-09-18_PSMfinalitems.csv")
df.dif.years.psm <- df.dif.years %>%
  select(any_of(final.psm.items))

df.dif.years.psm$F91 <- recode(df.dif.years.psm$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.dif.years.psm$F93 <- recode(df.dif.years.psm$F93, "4=3", as.factor = FALSE)
df.dif.years.psm$F98 <- recode(df.dif.years.psm$F98, "4=3", as.factor = FALSE)
df.dif.years.psm$F95 <- recode(df.dif.years.psm$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)

RIdifTable(df.dif.years.psm, dif.year)

Item 3 5 6 8 10 11 Mean location StDev MaxDiff
F88 -0.096 -0.077 -0.070 -0.043 -0.035 -0.045 -0.061 0.024 0.061
F91 0.104 0.149 0.102 0.132 0.078 0.060 0.104 0.033 0.088
F93 -0.144 -0.142 -0.125 -0.185 -0.211 -0.178 -0.164 0.032 0.085
F95 0.049 0.001 0.004 -0.009 0.038 0.002 0.014 0.024 0.058
F98 0.087 0.069 0.090 0.105 0.130 0.161 0.107 0.033 0.092
Code
RIdifFigure(df.dif.years.psm, dif.year)

Code
RIdifFigTime(df.dif.years.psm, dif.year)

Items är stabila och jämförbara över tid.

4.8 Reliabilitet - psykosomatiska besvär

Code
RItif(df.psm)

Vi har acceptabel reliabilitet för ca 68% av respondenterna.

4.9 Item-parametrar

Code
RIitemparams(df.dif.years.psm, "PsyksomItemParams.csv")
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Item location
F88 -1.09 -0.28 0.85 1.23 0.18
F91 0.01 0.70 NA NA 0.35
F93 -0.37 0.09 0.51 NA 0.08
F95 -0.05 0.22 0.57 NA 0.25
F98 -0.08 0.25 0.86 NA 0.34

4.10 Transformeringstabell

Med mätosäkerheter angivna som Logit std.error, där värden under 0.54 motsvarar reliabilitet över PSI = 0.7.

Code
RIscoreSE(df.dif.years.psm)
Ordinal sum score Logit score Logit std.error
0 -2.80 NA
1 -2.02 1.00
2 -1.32 0.72
3 -0.90 0.60
4 -0.58 0.53
5 -0.32 0.49
6 -0.09 0.47
7 0.13 0.46
8 0.34 0.46
9 0.56 0.47
10 0.79 0.49
11 1.04 0.53
12 1.35 0.59
13 1.77 0.71
14 2.45 0.99
15 3.20 NA

4.11 Rasch-analys 1 - psykiska besvär

Code
items.psm <- c("F88","F91","F93","F95","F98")

df.psy <- df.omit.na %>% 
  select(!any_of(items.psm))
itemnr item
F89 Känner du dig ledsen och deppig utan att veta varför?
F90 Händer det att du känner dig rädd utan att veta varför?
F92 Hur mycket skulle du vilja ändra på dig själv?
F94 Hur ofta tycker du att du inget duger till?
F96 Är du nöjd med ditt utseende?
F97 Känner du dig slö och olustig?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIitemfitPCM2(df.psy, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F89 0.773 0.79 -2.851 -2.865
F90 0.958 1.007 -0.342 -0.072
F92 0.778 0.784 -2.868 -3.095
F94 0.765 0.763 -2.414 -3.208
F96 0.879 0.882 -1.293 -1.324
F97 1.031 1.021 0.336 0.253
F99 1.059 1.016 0.664 0.155
Code
RIpcmPCA(df.psy)
PCA of Rasch model residuals
Eigenvalues
1.60
1.28
1.18
1.09
0.98
Code
RIresidcorr(df.psy, cutoff = 0.2)
F89 F90 F92 F94 F96 F97 F99
F89
F90 0.06
F92 -0.17 -0.18
F94 -0.11 -0.15 -0.08
F96 -0.34 -0.26 0.12 -0.05
F97 -0.14 -0.12 -0.24 -0.16 -0.2
F99 -0.12 -0.18 -0.18 -0.22 -0.11 -0.16
Note:
Relative cut-off value (highlighted in red) is 0.058, which is 0.2 above the average correlation.
Code
RItargeting(df.psy)

Code
RIitemHierarchy(df.psy)

Code
plot(mirt.rasch, type="trace")

Code
for (i in 1:ncol(df.psy)) {
    barplot(table(df.psy[, i]), col = "#8dc8c7", main = names(df.psy[i]), 
      ylab = "Number of responses", xlab = (itemlabels %>% 
                                              filter(itemnr %in% names(df.psy)) 
                                            %>% .[i,2]))
  }

Code
RIlistitems(df.psy)
itemnr item
F89 Känner du dig ledsen och deppig utan att veta varför?
F90 Händer det att du känner dig rädd utan att veta varför?
F92 Hur mycket skulle du vilja ändra på dig själv?
F94 Hur ofta tycker du att du inget duger till?
F96 Är du nöjd med ditt utseende?
F97 Känner du dig slö och olustig?
F99 Hur ofta tycker du att det är riktigt härligt att leva?

Många items har problem med oordnade svarskategorier.

Items 92 och 96 har en för stor residualkorrelation, men vi avvaktar med åtgärd tills svarskategorierna är åtgärdade.

4.12 Omkodning av svarskategorier

Vi behöver åtgärda flera items:

  • F89, 90, 96: 1+2 och 3+4
  • F92, 99: 3+4
  • F97: 0+1

(flera varianter testades, och ovanstående fungerade minst dåligt)

Code
df.psy$F89 <- recode(df.psy$F89, "2=1;3=2;4=2", as.factor = FALSE)
df.psy$F90 <- recode(df.psy$F90, "2=1;3=2;4=2", as.factor = FALSE)
df.psy$F96 <- recode(df.psy$F96, "2=1;3=2;4=2", as.factor = FALSE)
df.psy$F97 <- recode(df.psy$F97, "1=0;2=1;3=2;4=3", as.factor = FALSE)

df.psy$F92 <- recode(df.psy$F92, "4=3", as.factor = FALSE)
df.psy$F99 <- recode(df.psy$F99, "4=3", as.factor = FALSE)
Code
plot(mirt.rasch, type="trace")

Code
RIitemHierarchy(df.psy)

4.13 Rasch-analys 2 - psykiska besvär

Efter att ha åtgärdat svarskategorierna.

itemnr item
F89 Känner du dig ledsen och deppig utan att veta varför?
F90 Händer det att du känner dig rädd utan att veta varför?
F92 Hur mycket skulle du vilja ändra på dig själv?
F94 Hur ofta tycker du att du inget duger till?
F96 Är du nöjd med ditt utseende?
F97 Känner du dig slö och olustig?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIitemfitPCM2(df.psy, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F89 0.777 0.79 -3.002 -3.035
F90 0.912 0.951 -0.68 -0.477
F92 0.832 0.841 -1.909 -2.269
F94 0.774 0.792 -2.438 -2.786
F96 0.845 0.845 -2.073 -2.126
F97 1.071 1.024 0.707 0.404
F99 1.076 1.045 0.762 0.487
Code
RIpcmPCA(df.psy)
PCA of Rasch model residuals
Eigenvalues
1.52
1.33
1.21
1.13
0.91
Code
RIresidcorr(df.psy, cutoff = 0.2)
F89 F90 F92 F94 F96 F97 F99
F89
F90 0.12
F92 -0.13 -0.13
F94 -0.11 -0.13 -0.17
F96 -0.2 -0.18 0.12 -0.07
F97 -0.11 -0.11 -0.24 -0.22 -0.17
F99 -0.11 -0.16 -0.24 -0.31 -0.12 -0.16
Note:
Relative cut-off value (highlighted in red) is 0.065, which is 0.2 above the average correlation.
Code
RItargeting(df.psy)

Code
RIitemHierarchy(df.psy)

Vi har två residualkorrelationer ca 0.25 över medelvärdet för samtliga residualkorrelationer.

  • F96 tas bort (F92 behålls pga bättre targeting)
  • F90 tas bort (F89 behålls pga bättre targeting)

4.13.1 Reliabilitet

Även om vi inte undersökt DIF, vilket bör göras innan reliabiliteten bedöms, så har de påtagliga problemen med svarskategorierna medfört att vi kan ha såpass låg reliabilitet att det inte är meningsfullt med DIF-analysen.

Code
df.psy$F96 <- NULL
df.psy$F90 <- NULL

RItif(df.psy)

Eftersom många items uppvisade problem med svarskategorierna blir reliabiliteten låg.

Vi provar att slå samman kvarvarande items från psykiska besvär med de items som utgör psykosomatiska besvär.

Code
df.psf <- cbind(df.psm,df.psy)

df.psf <- df.psf |> 
  relocate(F92, .after = "F91") |> 
  relocate(F97, .before = "F98") |> 
  relocate(F89, .after = "F88") |>
  relocate(F94, .before = "F95")

# psf.items <- names(df.psf)
# 
# df.psf <- df.omit.na |> 
#   select(any_of(psf.items))

4.14 Rasch-analys 1 psykiska/psykosomatiska besvär

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F89 Känner du dig ledsen och deppig utan att veta varför?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94 Hur ofta tycker du att du inget duger till?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIitemfitPCM2(df.psf, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F88 1.078 1.058 1.263 0.821
F89 0.719 0.728 -4.626 -4.358
F91 0.963 0.961 -0.507 -0.437
F92 0.86 0.863 -1.722 -2.189
F93 1.014 0.993 0.468 -0.336
F94 0.904 0.887 -1.003 -1.456
F95 1.006 0.985 -0.279 0.111
F97 0.967 0.954 -0.201 -0.423
F98 0.95 0.957 -0.996 -0.749
F99 1.064 1.01 0.696 -0.026
Code
RIpcmPCA(df.psf)
PCA of Rasch model residuals
Eigenvalues
1.78
1.26
1.13
1.08
1.02
Code
RIloadLoc(df.psf)

Code
RIresidcorr(df.psf, cutoff = 0.2)
F88 F89 F91 F92 F93 F94 F95 F97 F98 F99
F88
F89 -0.1
F91 -0.07 -0.03
F92 -0.22 0.05 -0.06
F93 -0.03 -0.07 0 -0.1
F94 -0.23 0.08 -0.16 0.11 -0.18
F95 -0.08 -0.14 -0.04 -0.21 -0.11 -0.26
F97 -0.14 0.01 -0.13 -0.05 -0.18 0.02 -0.13
F98 -0.05 -0.09 -0.08 -0.2 -0.1 -0.21 0.11 -0.17
F99 -0.23 0.04 -0.12 0.01 -0.2 -0.02 -0.14 0.01 -0.19
Note:
Relative cut-off value (highlighted in red) is 0.109, which is 0.2 above the average correlation.
Code
RItargeting(df.psf)

Code
RIitemHierarchy(df.psf)

Code
plot(mirt.rasch, type="trace")

Code
df.psf$F89 <- NULL
df.psf$F94 <- recode(df.psf$F94, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.erm <- PCM(df.psf)
plotICC(df.erm, item.subset = "F94")

Code
RIitemHierarchy(df.psf)

F94 har fortfarande problem med svarskategorierna. Vi slår samman de två mittersta.

Code
df.psf$F94 <- recode(df.psf$F94, "3=2;2=1", as.factor = FALSE)
df.erm <- PCM(df.psf)
plotICC(df.erm, item.subset = "F94")

Code
RIitemHierarchy(df.psf)

4.15 Rasch-analys 2 psykiska/psykosomatiska besvär

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94 Hur ofta tycker du att du inget duger till?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIitemfitPCM2(df.psf, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F88 1.009 0.997 0.038 -0.06
F91 0.929 0.936 -0.919 -1.045
F92 0.871 0.874 -1.736 -1.727
F93 0.961 0.95 -0.4 -0.571
F94 0.835 0.839 -1.932 -2.329
F95 0.922 0.918 -0.992 -1.229
F97 0.959 0.948 -0.314 -0.524
F98 0.895 0.908 -0.992 -1.489
F99 1.058 1.009 0.972 0.213
Code
RIpcmPCA(df.psf)
PCA of Rasch model residuals
Eigenvalues
1.69
1.29
1.16
1.07
1.01
Code
RIloadLoc(df.psf)

Code
RIresidcorr(df.psf, cutoff = 0.2)
F88 F91 F92 F93 F94 F95 F97 F98 F99
F88
F91 -0.09
F92 -0.22 -0.05
F93 -0.06 -0.02 -0.1
F94 -0.14 -0.09 0.16 -0.11
F95 -0.12 -0.07 -0.21 -0.15 -0.16
F97 -0.15 -0.13 -0.03 -0.18 0.09 -0.15
F98 -0.08 -0.1 -0.2 -0.13 -0.12 0.07 -0.18
F99 -0.24 -0.12 0.03 -0.2 0.05 -0.16 0.01 -0.2
Note:
Relative cut-off value (highlighted in red) is 0.102, which is 0.2 above the average correlation.
Code
RItargeting(df.psf)

Code
RIitemHierarchy(df.psf)

Code
plot(mirt.rasch, type="trace")

F94 korrelerar för starkt med F92. Vi tar bort F94 p.g.a. sämre targeting och ZSTD.

Code
df.psf$F94 <- NULL

4.16 Rasch-analys 3 psykiska/psykosomatiska besvär

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIitemfitPCM2(df.psf, 300, 32, 8)
OutfitMSQ InfitMSQ OutfitZSTD InfitZSTD
F88 0.962 0.958 -0.554 -0.566
F91 0.902 0.915 -1.238 -1.462
F92 0.886 0.885 -1.667 -1.734
F93 0.923 0.922 -0.612 -1.136
F95 0.882 0.881 -1.355 -1.819
F97 0.972 0.947 -0.68 -0.84
F98 0.863 0.878 -1.736 -1.777
F99 1.055 1.004 0.646 0.296
Code
RIpcmPCA(df.psf)
PCA of Rasch model residuals
Eigenvalues
1.55
1.31
1.17
1.05
1.00
Code
RIloadLoc(df.psf)

Code
RIresidcorr(df.psf, cutoff = 0.2)
F88 F91 F92 F93 F95 F97 F98 F99
F88
F91 -0.1
F92 -0.21 -0.04
F93 -0.07 -0.02 -0.08
F95 -0.14 -0.08 -0.2 -0.16
F97 -0.14 -0.13 0 -0.17 -0.14
F98 -0.09 -0.11 -0.18 -0.14 0.06 -0.17
F99 -0.24 -0.12 0.04 -0.2 -0.16 0.03 -0.2
Note:
Relative cut-off value (highlighted in red) is 0.087, which is 0.2 above the average correlation.
Code
RItargeting(df.psf)

Code
RIitemHierarchy(df.psf)

Detta ser ut att fungera bra. Vi övergår till DIF-analys

4.17 DIF/invarians

4.17.1 Kön

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
#dif.gender<-recode(dif.gender,"'Pojke'=1;'Flicka'=2",as.factor=FALSE)
RIdifTable(df.psf, dif.gender)

Item 2 3 Mean location StDev MaxDiff
F88 -0.204 -0.012 -0.108 0.136 0.193
F91 0.056 0.086 0.071 0.021 0.030
F92 -0.623 -0.443 -0.533 0.127 0.180
F93 -0.308 -0.180 -0.244 0.091 0.128
F95 -0.002 -0.171 -0.086 0.119 0.169
F97 0.779 0.628 0.704 0.106 0.151
F98 -0.048 0.056 0.004 0.074 0.104
F99 0.350 0.035 0.193 0.223 0.315
Code
RIdifFigure(df.psf, dif.gender)

Code
RIdifFigThresh(df.psf, dif.gender)

F99 är enda item som går strax över 0.3 logits. Behöver ej åtgärdas.

4.17.2 Årskurs

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIdifTable(df.psf, dif.arskurs)

Item 2 3 Mean location StDev MaxDiff
F88 -0.157 -0.077 -0.117 0.056 0.079
F91 0.004 0.111 0.057 0.076 0.107
F92 -0.580 -0.517 -0.548 0.044 0.062
F93 -0.208 -0.286 -0.247 0.055 0.078
F95 -0.090 -0.085 -0.088 0.003 0.005
F97 0.744 0.699 0.722 0.032 0.045
F98 0.086 -0.058 0.014 0.102 0.144
F99 0.200 0.214 0.207 0.010 0.014
Code
RIdifFigure(df.psf, dif.arskurs)

Code
RIdifFigThresh(df.psf, dif.arskurs)

Inga skillnader mellan årskurser.

4.17.3 Årtal

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
final.items.psf <- names(df.psf)
#write.csv(final.items, file = "2022-09-16 IFoptimalItems.csv")
df.dif.years <- df.dif.years %>% 
  select(any_of(final.items.psf))

df.dif.years$F91 <- recode(df.dif.years$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.dif.years$F92 <- recode(df.dif.years$F92, "4=3", as.factor = FALSE)
df.dif.years$F93 <- recode(df.dif.years$F93, "4=3", as.factor = FALSE)
df.dif.years$F98 <- recode(df.dif.years$F98, "4=3", as.factor = FALSE)
df.dif.years$F99 <- recode(df.dif.years$F99, "4=3", as.factor = FALSE)
df.dif.years$F95 <- recode(df.dif.years$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.dif.years$F97 <- recode(df.dif.years$F97, "4=3;3=2;2=1;1=0", as.factor = FALSE)

RIdifTable(df.dif.years, dif.year)

Item 3 6 7 8 11 12 14 15 Mean location StDev MaxDiff
F88 -0.179 -0.177 -0.180 -0.161 -0.111 -0.148 -0.112 -0.141 -0.151 0.028 0.070
F91 0.017 0.046 0.039 0.007 0.060 0.022 -0.003 -0.036 0.019 0.030 0.096
F92 -0.502 -0.504 -0.494 -0.518 -0.542 -0.508 -0.495 -0.503 -0.508 0.016 0.048
F93 -0.221 -0.241 -0.233 -0.210 -0.247 -0.277 -0.278 -0.266 -0.247 0.025 0.068
F95 -0.038 -0.121 -0.082 -0.088 -0.088 -0.102 -0.042 -0.092 -0.082 0.028 0.084
F97 0.777 0.811 0.781 0.744 0.719 0.769 0.693 0.741 0.754 0.038 0.118
F98 -0.003 -0.045 -0.029 -0.007 0.002 0.021 0.045 0.061 0.006 0.035 0.106
F99 0.149 0.232 0.198 0.234 0.208 0.223 0.193 0.234 0.209 0.029 0.086
Code
RIdifFigure(df.dif.years, dif.year)

Code
RIdifFigTime(df.dif.years, dif.year)

Code
RIdifFigThresh(df.dif.years, dif.year)

Items fungerar stabilt över tid.

4.17.4 Interaktion årskurs och kön

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
dfin <- df.psf

df.tree <- data.frame(matrix(ncol = 0, nrow = nrow(dfin)))
df.tree$difdata <- as.matrix(dfin)
df.tree$dif.gender <- dif.gender
df.tree$dif.arskurs <- dif.arskurs

pctree.out <- pctree(difdata ~ dif.gender + dif.arskurs, data = df.tree)
plot(pctree.out)

Code
cutoff <- 0.5
diffig <- itempar(pctree.out) %>%
  as.data.frame() %>%
  t() %>%
  as.data.frame() %>%
  mutate(
    `Mean location` = rowMeans(.),
    StDev = rowSds(as.matrix(.))
  ) %>%
  rowwise() %>%
  mutate(MaxDiff = (max(c_across(c(1:(ncol(.) - 2))))) -
    min(c_across(c(1:(ncol(.) - 2))))) %>%
  ungroup() %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  rownames_to_column(var = "Item") %>%
  mutate(Item = names(dfin)) %>%
  relocate(MaxDiff, .after = last_col())

  formattable(diffig,
    list(MaxDiff = formatter("span",
      style = ~ style(color = ifelse(MaxDiff < -cutoff,
        "red", ifelse(MaxDiff > cutoff, "red", "black")
      ))
    ),
    formattable::area(col = 2:6) ~ color_tile(RISEprimYellowLight, RISEprimGreen)
    ),
    table.attr = "class=\"table table-striped\" style=\"font-size: 15px; font-family: Lato\""
  )
Item 3 4 6 7 Mean location StDev MaxDiff
F88 -0.235 -0.184 -0.071 0.036 -0.113 0.121 0.271
F91 -0.029 0.114 0.046 0.119 0.062 0.070 0.149
F92 -0.692 -0.576 -0.439 -0.446 -0.538 0.120 0.253
F93 -0.241 -0.355 -0.167 -0.190 -0.238 0.083 0.188
F95 0.019 -0.018 -0.183 -0.160 -0.085 0.101 0.202
F97 0.820 0.753 0.627 0.628 0.707 0.096 0.193
F98 0.075 -0.133 0.092 0.028 0.016 0.103 0.225
F99 0.283 0.399 0.095 -0.015 0.190 0.186 0.414

4.17.5 Item fit över tid

DIF-analyserna visar på stabilitet i item locations över tid. Det kan också vara av intresse att se hur item fit ser ut över tid (2006-2020).

Code
final.items <- read_csv("PSFitemnr.csv") %>% 
  pull(itemnr)

df.fit <- df.all %>%
  select(final.items,ar) %>% 
  na.omit()

df.fit$F91 <- recode(df.fit$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.fit$F92 <- recode(df.fit$F92, "4=3", as.factor = FALSE)
df.fit$F93 <- recode(df.fit$F93, "4=3", as.factor = FALSE)
df.fit$F98 <- recode(df.fit$F98, "4=3", as.factor = FALSE)
df.fit$F99 <- recode(df.fit$F99, "4=3", as.factor = FALSE)
df.fit$F95 <- recode(df.fit$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.fit$F97 <- recode(df.fit$F97, "4=3;3=2;2=1;1=0", as.factor = FALSE)
Code
library(furrr)
plan(multisession, workers = 8)

itemfit.time <- df.fit %>% 
  split(.$ar) %>% 
  future_map(~ RIitemfitPCM(subset(., select = -ar), table = FALSE), 
             .options = furrr_options(seed = TRUE))

årtal <- c("2006", "2008", "2010", "2012", "2014", "2016", "2018", "2020")
itemfit.all <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(itemfit.all) <- c("Item", "År", 'InfitMSQ','OutfitMSQ')

for (i in årtal) {
  itemfit.add <- itemfit.time[[i]] %>% 
    rownames_to_column("Item") %>% 
    add_columnr = i)
  
  itemfit.all <- rbind(itemfit.all, itemfit.add)
}
Code
itemfit.all %>%
  group_by(Item) %>%
  summarise(
    mean_Out = mean(OutfitMSQ),
    sd_Out = sd(OutfitMSQ),
    max_Out = max(OutfitMSQ),
    min_Out = min(OutfitMSQ),
    mean_In = mean(InfitMSQ),
    sd_In = sd(InfitMSQ),
    max_In = max(InfitMSQ),
    min_In = min(InfitMSQ)
  ) %>%
  mutate(
    diffmax_OutfitMSQ = max_Out - min_Out,
    diffmax_InfitMSQ = max_In - min_In,
    item = Item
  ) %>%
  mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
  mutate(across(starts_with(c("max", "min")), ~ cell_spec(.x, color = ifelse(.x < msq_min, "red",
    ifelse(.x > msq_max, "red", "black")
  )))) %>%
  kbl(booktabs = T, escape = F,
      col.names = c("Item","Mean","SD","Max","Min","Mean","SD","Max","Min","OutfitMSQ","InfitMSQ","Item")) %>%
  # bootstrap options are for HTML output
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    position = "left",
    full_width = F,
    font_size = 15,
    fixed_thead = T
  ) %>%
  column_spec(c(1, 12), bold = T) %>%
  kable_classic(html_font = "Lato") %>%
  # latex_options are for PDF output
  kable_styling(latex_options = c("striped", "scale_down")) %>%
  add_header_above(c(" " = 1, 
                     "OutfitMSQ" = 4, 
                     "InfitMSQ" = 4, 
                     "Max differences" = 2,
                     " " = 1))
OutfitMSQ
InfitMSQ
Max differences
Item Mean SD Max Min Mean SD Max Min OutfitMSQ InfitMSQ Item
F88 0.97 0.02 1 0.96 0.96 0.01 0.99 0.95 0.04 0.04 F88
F91 0.88 0.01 0.91 0.86 0.89 0.01 0.92 0.87 0.05 0.04 F91
F92 0.91 0.02 0.95 0.89 0.91 0.02 0.94 0.89 0.06 0.05 F92
F93 0.93 0.01 0.95 0.91 0.93 0.01 0.94 0.91 0.04 0.03 F93
F95 0.88 0.03 0.94 0.85 0.89 0.02 0.93 0.86 0.09 0.07 F95
F97 0.95 0.03 0.99 0.9 0.94 0.02 0.96 0.9 0.09 0.06 F97
F98 0.86 0.01 0.86 0.84 0.87 0.01 0.88 0.86 0.02 0.02 F98
F99 1.04 0.02 1.06 1.01 1.00 0.01 1.01 0.98 0.05 0.03 F99

Inget item går utanför gränsvärden 0.7-1.3.

4.18 Reliabilitet - psykiska/psykosomatiska besvär

Code
RItif(df.psf)

Bättre reliabilitet än vad någon av delskalorna uppnår separat.

4.19 Person fit

Code
RIpfit(df.psf)

4.20 Item-parametrar

itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
RIitemparams(df.dif.years)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Item location
F88 -1.09 -0.37 0.74 1.01 0.07
F91 -0.06 0.55 NA NA 0.25
F92 -1.36 0.06 0.46 NA -0.28
F93 -0.41 -0.01 0.36 NA -0.02
F95 -0.10 0.13 0.41 NA 0.14
F97 0.36 1.06 1.52 NA 0.98
F98 -0.14 0.14 0.68 NA 0.23
F99 -0.32 0.60 1.03 NA 0.44
Code
itemlabels %>% 
  filter(itemnr %in% names(df.dif.years)) %>% 
  write_csv("PSFitemnr.csv")
Code
RIscoreSE(df.dif.years)
Ordinal sum score Logit score Logit std.error
0 -3.38 NA
1 -2.59 1.01
2 -1.88 0.72
3 -1.45 0.60
4 -1.14 0.52
5 -0.89 0.47
6 -0.68 0.44
7 -0.50 0.41
8 -0.34 0.40
9 -0.19 0.38
10 -0.04 0.38
11 0.10 0.37
12 0.24 0.37
13 0.37 0.37
14 0.51 0.38
15 0.66 0.38
16 0.81 0.39
17 0.97 0.41
18 1.14 0.43
19 1.34 0.46
20 1.58 0.51
21 1.87 0.58
22 2.27 0.70
23 2.94 0.99
24 3.69 NA

4.21 Exempel med single-item vs index

Data från enbart 2020.

OBS! Viktigt att komma ihåg att höga värden = hög risk, vilket även gäller svar på enskilda frågor. En positiv fråga som F99 har alltså omvända svarskategorier (3 och 4 sammanslagna):

  • 3-4 - Sällan / Någon enstaka gång
  • 2 - Ibland
  • 1 - Ganska ofta
  • 0 - Väldigt ofta
Code
#---- create df with only 2020 respondents as example----
df.if <- df %>% 
  filter(ar == 2020) %>% 
  select(any_of(names(df.psf)),Kön) %>% 
  na.omit()

F99o <- df.if$F99

# rcat2 <- c("F95","F97")
# rcat3 <- c("F92","F93","F99")
# 
# for (i in rcat2) {
#   df.if[[i]]<-recode(df.if[[i]],"1=0;2=1;3=2;4=3",as.factor=FALSE)
# }
# for (i in rcat3) {
#   df.if[[i]]<-recode(df.if[[i]],"4=3",as.factor=FALSE)
# }
df.if$F91 <- recode(df.if$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.if$F92 <- recode(df.if$F92, "4=3", as.factor = FALSE)
df.if$F93 <- recode(df.if$F93, "4=3", as.factor = FALSE)
df.if$F95 <- recode(df.if$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.if$F97 <- recode(df.if$F97, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.if$F98 <- recode(df.if$F98, "4=3", as.factor = FALSE)
df.if$F99 <- recode(df.if$F99, "4=3", as.factor = FALSE)
df.if.gender <- df.if$Kön
df.if$Kön <- NULL
#---- estimate person locations----
library(catR)

# create matrix with item parameters
df.erm <- PCM(df.psf)
item.estimates <- eRm::thresholds(df.erm)
item_difficulty <- as.data.frame(item.estimates[["threshtable"]][["1"]])
item_difficulty$Location <- NULL
items <- item_difficulty %>% 
  mutate_if(is.character, as.numeric) %>% 
  as.matrix()

# loop through all participants to calculate theta for each one
thetaEstScores <- c()
#p1 <- c()
for (i in 1:nrow(df.if)){
  p1 <- as.numeric(as.vector(df.if[i,]))
  ptheta <- thetaEst(items, p1, model = "PCM", method = "ML") # several method options available for theta calculation
  thetaEstScores <- c(thetaEstScores,ptheta)
}

# # loop through all participants theta scores to calculate SE for each one
# thetaEstSE <- c()
# for (i in 1:length(thetaEstScores)) {
#   p1theta <- thetaEstScores[i]
#   p1 <- as.numeric(as.vector(df.if[i,]))
#   pthetaSE <- semTheta(thEst = p1theta, it = items, x = p1, model = "PCM", method = "ML") # several method options available for theta calculation
#   thetaEstSE <- c(thetaEstSE,pthetaSE)
# }

df.if$SumScores <- thetaEstScores
df.if$Kön <- df.if.gender
df.if$F99o <- F99o
#write_parquet(df.if, sink = "PSF2020.parquet")

#---- barplot for single item----
barplot(table(df.if$F99), col = "#8dc8c7", main = "Svarsfördelning för item F99")

Code
#---- visualize single item vs index score----
# some code from https://dallasnova.rbind.io/post/efficient-data-visualization-with-faded-raincloud-plots-delete-boxplot/
library(ggdist) # for shadeable density slabs
library(gghalves) # for half-half geoms
library(ggpp) # for position_dodge2nudge
library(colorspace) # for lightening color palettes
library(extrafont) # for Lato, RISE font
library(stringr)
# Setting colorblind-friendly palette
cbPalette <-c("#999999","#E69F00", "#56B4E9","#009E73",
              "#F0E442", "#0072B2", "#D55E00","#CC79A7")

cbPalette <- lighten(cbPalette, amount = 0.6, space = "HLS")

#---- plot for index value distributions divided by gender----

ggplot(df.if, aes(x = Kön, y = SumScores)) +
  stat_slab(side = "right", show.legend = F,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=Kön), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal",show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  #xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Svarsfördelning för psykisk och psykosomatisk hälsa") +
  #labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
  #     subtitle = "Item F99 - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

Code
#---- plots for index value vs single item response ----
ggplot(df.if, aes(x = as.factor(F99), y = SumScores)) +
  #geom_violin()
  # density distribution slab
  #geom_boxplot(size = 0.2) +
  stat_slab(side = "right", show.legend = F,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=as.factor(F99)), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal", show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Jämförelse av ett item och indexvärde") +
  labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
       subtitle = "Item F99 - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

Code
ggplot(df.if, aes(x = as.factor(F99), y = SumScores)) +
  #geom_violin()
  # density distribution slab
  #geom_boxplot(size = 0.2) +
  stat_slab(side = "right", show.legend = T,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=Kön), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal",show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Jämförelse av ett item och indexvärde") +
  labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
       subtitle = "Item F99 - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

Code
# stat_dotsinterval(scale = 0.5, quantiles = 100, position = "dodge") +

4.21.1 Distribution med ej omkodade svarskategorier

Code
barplot(table(df.if$F99o), col = "#8dc8c7", main = "Svarsfördelning för item F99 innan recode")

Code
#plotICC(df.ermo, item.subset = "F99")

ggplot(df.if, aes(x = as.factor(F99o), y = SumScores)) +
  #geom_violin()
  # density distribution slab
  #geom_boxplot(size = 0.2) +
  stat_slab(side = "right", show.legend = F,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=as.factor(F99o)), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal", show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Jämförelse av ett item och indexvärde") +
  labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
       subtitle = "Item F99o - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

4.22 Jämförelser befintligt och nytt index

Rapporter utifrån Stockholmsenkäten har hittills räknat samman ett index för “Psykisk hälsa”, baserat på enkätsvar för de fem items som i denna analys återfinns under rubriken “Psykosomatiska besvär”

Code
items.psm <- c("F88","F91","F93","F95","F98")

df.index <- df.all %>%
  select(items.psm, Kön, ARSKURS, ar, SkolSDO)

RIlistitems(df.index)
itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?

Indexberäkningen som använts utgår från att svaren översätts till siffror och sedan beräknas ett medelvärde av alla items. Lägsta svarsalternativet (“Aldrig”) ersätts med siffran 100, näst lägsta med 75, och så vidare. Det innebär att ett högt värde motsvarar låg förekomst av besvär.

I Stockholm Stads rapport för år 2022 återges resultatet på detta vis:

Vi kan återskapa figuren utifrån data från 2020 för årskurs 9.

Code
df.indexR <- df.index %>%
  mutate(across(items.psm, ~ recode(.x, "0=100;1=75;2=50;3=25;4=0"))) %>%
  na.omit() %>%
  mutate(PSFindex = rowMeans(across(items.psm)))

library(ggchicklet)
# create a wide palette based on RISE colors
RISEpalette1 <- colorRampPalette(colors = c("#009ca6", "#e83c63", "#ffe500"))(6)
# "#009CA6" "#5C758B" "#B94F70" "#EC5D4F" "#F5A127" "#FFE500"
RISEpalette2 <- colorRampPalette(colors = c("#009ca6", "#e83c63", "#ffe500"))(8)

df.indexR %>%
  filter(ar == 2020) %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PSFindex, na.rm = T), 0)) %>%
  ggplot(aes(y = PSFmean, x = SkolSDO, fill = Kön, color = Kön)) +
  geom_chicklet(
    radius = grid::unit(1.5, "mm"),
    position = "dodge"
  ) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), linetype = 2, color = "darkgrey") +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), linetype = 3, color = "black") +
  geom_text(aes(label = PSFmean, y = PSFmean + 2),
    position = position_dodge(width = .9),
    color = "black", size = 3.8
  ) +

  # annotate("text", label = "Medelvärde", y = 57, x = 3, color = "darkgrey") +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_discrete("Stadsdel", guide = guide_axis(n.dodge = 3)) +
  scale_y_continuous(limits = c(0,100)) +
  theme_minimal(base_size = 12) +
  theme_rise(axissize = 14, titlesize = 16) +
  theme(text = element_text(family = "Lato")) +
  ylab("Medelvärde") +
  labs(
    title = "Psykisk hälsa indexvärde",
    subtitle = "Befintligt index, 2020, åk 9",
    caption = "Streckade linjer betecknar medelvärde för samtliga stadsdelar"
  )

4.22.1 Nytt index (5 items)

Här utgår vi från samma fem items som i det befintliga indexet, men kodar om svarskategorierna som inte fungerade som tänkt, och estimerar fram ett indexvärde för varje individ utifrån Rasch-parametrarna. För jämförbarhet skalas sedan mätvärdet om till att omfånget 0-100.

Frågorna mäter i huvudsak besvär, vilket medfört att indexet uppkallats efter det. För jämförbarhet med det befintliga indexet “psykisk hälsa”, har mätvärden reverserats, så högre värden = lägre nivåer av besvär. Samma gäller indexet med 8 items.

Code
# read item parameters from previous analysis
itemParams.psm <- read_csv("PsyksomItemParams.csv") %>% 
  as.matrix()

# create id variable to re-join thetas
df.index <- df.index %>% 
  mutate(id = seq.int(nrow(.)))

# recode items to correct response categories
df.indexNew <- df.index %>% 
  select(all_of(c(items.psm,"id"))) %>% 
  na.omit()

df.indexNew$F91 <- recode(df.indexNew$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.indexNew$F93 <- recode(df.indexNew$F93, "4=3", as.factor = FALSE)
df.indexNew$F98 <- recode(df.indexNew$F98, "4=3", as.factor = FALSE)
df.indexNew$F95 <- recode(df.indexNew$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)

#RItileplot(df.indexNew %>% select(!id))
#itemParams.psm
df.indexNew$PSMscore <- df.indexNew %>% 
  select(!id) %>% 
  RIestThetas2(cpu = 8, itemParams = itemParams.psm)
df.index <- df.indexNew %>% 
  select(id,PSMscore) %>% 
  left_join(df.index,df.indexNew, by = "id")

# reverse score to match old index
df.index$PSMscore <- df.index$PSMscore*-1
Code
# rescale to match 0-100 range
df.index$PSMscore100 <- round(scales::rescale(df.index$PSMscore, to = c(0,100)),0)

4.22.2 Nytt index (8 items)

Det index som både representerar psykosomatiska och psykiska besvär innehåller 8 frågor, och beräknas utifrån items efter korrigering av svarskategorier.

Code
newIndex <- read_csv("PSFitemnr.csv")

items.psnew <- newIndex$itemnr

df.index8 <- df.all %>%
  select(items.psnew, Kön, ARSKURS, ar, SkolSDO)

RIlistitems(df.index8)
itemnr item
F88 Hur ofta har du haft huvudvärk detta läsår?
F91 Hur ofta har du dålig aptit?
F92 Hur mycket skulle du vilja ändra på dig själv?
F93 Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95 Hur ofta har du under detta läsår haft svårt att somna?
F97 Känner du dig slö och olustig?
F98 Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99 Hur ofta tycker du att det är riktigt härligt att leva?
Code
# read item parameters from previous analysis
itemParams.psnew <- read_csv("itemParameters.csv") %>% 
  as.matrix()

# create id variable to re-join thetas
df.index8 <- df.index8 %>% 
  mutate(id = seq.int(nrow(.)))

# recode items to correct response categories
df.indexNew8 <- df.index8 %>% 
  select(all_of(c(items.psnew,"id"))) %>% 
  na.omit()

rcat1 <- c("F91")
rcat2 <- c("F95", "F97")
rcat3 <- c("F92", "F93", "F98", "F99")

for (i in rcat1) {
  df.indexNew8[[i]] <- recode(df.indexNew8[[i]], "2=1;3=2;4=2", as.factor = FALSE)
}
for (i in rcat2) {
  df.indexNew8[[i]] <- recode(df.indexNew8[[i]], "1=0;2=1;3=2;4=3", as.factor = FALSE)
}
for (i in rcat3) {
  df.indexNew8[[i]] <- recode(df.indexNew8[[i]], "4=3", as.factor = FALSE)
}

#RItileplot(df.indexNew %>% select(!id))
#itemParams.psm
df.indexNew8$PS8score <- df.indexNew8 %>% 
  select(!id) %>% 
  RIestThetas2(cpu = 8, itemParams = itemParams.psnew)
df.index8 <- df.indexNew8 %>% 
  select(id,PS8score) %>% 
  left_join(df.index8,df.indexNew8, by = "id")

# reverse score to match old index
df.index8$PS8score <- df.index8$PS8score*-1
# and rescaled version
df.index8$PS8score100 <- round(scales::rescale(df.index8$PS8score, to = c(0,100)),0)
Code
df.index %>% 
  filter(ar == 2020) %>% 
  filter(Kön %in% c("Pojke","Flicka")) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  filter(!is.na(SkolSDO)) %>% 
  group_by(SkolSDO,Kön) %>% 
  summarise(PSFmean = round(mean(PSMscore100, na.rm = T),0)) %>% 
  ggplot(aes(y = PSFmean, x = SkolSDO, fill = Kön, color = Kön)) +
  geom_chicklet(
    radius = grid::unit(1.5, "mm"),
    position = "dodge"
  ) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), linetype = 2, color = "darkgrey") +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), linetype = 3, color = "black") +
  geom_text(aes(label = PSFmean, y = PSFmean + 2),
    position = position_dodge(width = .9),
    color = "black", size = 3.8
  ) +

  # annotate("text", label = "Medelvärde", y = 57, x = 3, color = "darkgrey") +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_discrete("Stadsdel", guide = guide_axis(n.dodge = 3)) +
  scale_y_continuous(limits = c(0,100)) +
  theme_minimal(base_size = 12) +
  theme_rise(axissize = 14, titlesize = 16) +
  theme(text = element_text(family = "Lato")) +
  ylab("Medelvärde") +
  labs(title = "Psykosomatiska besvär indexvärde",
       subtitle = "5 items, år 2020, åk 9",
       caption = "Streckade linjer betecknar medelvärde per kön för samtliga stadsdelar")

Code
df.index8 %>% 
  filter(ar == 2020) %>% 
  filter(Kön %in% c("Pojke","Flicka")) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  filter(!is.na(SkolSDO)) %>% 
  group_by(SkolSDO,Kön) %>% 
  summarise(PSFmean = round(mean(PS8score100, na.rm = T),0)) %>% 
  ggplot(aes(y = PSFmean, x = SkolSDO, fill = Kön, color = Kön)) +
  geom_chicklet(
    radius = grid::unit(1.5, "mm"),
    position = "dodge"
  ) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), linetype = 2, color = "darkgrey") +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), linetype = 3, color = "black") +
  geom_text(aes(label = PSFmean, y = PSFmean + 2),
    position = position_dodge(width = .9),
    color = "black", size = 3.8
  ) +

  # annotate("text", label = "Medelvärde", y = 57, x = 3, color = "darkgrey") +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_discrete("Stadsdel", guide = guide_axis(n.dodge = 3)) +
  scale_y_continuous(limits = c(0,100)) +
  theme_minimal(base_size = 12) +
  theme_rise(axissize = 14, titlesize = 16) +
  theme(text = element_text(family = "Lato")) +
  ylab("Medelvärde") +
  labs(title = "Psykiska/psykosomatiska besvär indexvärde",
       subtitle = "8 items, år 2020, åk 9",
       caption = "Streckade linjer betecknar medelvärde per kön för samtliga stadsdelar")

4.22.3 Indexvärde över tid, åk 9

Code
årtal <- c(2006,2008,2010,2012,2014,2016,2018,2020,2022)
# library(ggiraph)
# psfTid <- df.indexR %>%
#   filter(Kön %in% c("Pojke", "Flicka")) %>%
#   filter(!is.na(SkolSDO)) %>%
#   group_by(ar, SkolSDO, Kön) %>%
#   summarise(PSFmean = round(mean(PSFindex, na.rm = T), 0),
#             stdev = round(sd(PSFindex, na.rm = T))) %>%
#   ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
#   geom_line(linewidth = 1) +
#   geom_point_interactive(aes(tooltip = PSFmean),
#                          size = 1.5) +
#   geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
#               alpha = 0.15, linetype = 0) +
#   geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
#              linetype = 2, color = "black", alpha = 0.8) +
#   geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
#              linetype = 3, color = "black", alpha = 0.8) +
#   scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
#   scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
#   scale_y_continuous("Medelvärde", limits = c(0,100)) +
#   labs(title = "Psykisk ohälsa indexvärde",
#        caption = "Skuggade fält indikerar en standardavvikelse") +
#   facet_wrap(~SkolSDO) +
#   theme_minimal(base_size = 11) +
#   theme_rise()
# 
# ggiraph(ggobj = psfTid)

df.indexR %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(ar, SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PSFindex, na.rm = T), 0),
            stdev = round(sd(PSFindex, na.rm = T))) %>%
  ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
              alpha = 0.15, linetype = 0) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
             linetype = 2, color = "black", alpha = 0.8) +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
             linetype = 3, color = "black", alpha = 0.8) +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
  scale_y_continuous("Medelvärde", limits = c(0,100)) +
  labs(title = "Psykisk hälsa indexvärde",
        subtitle = "5 items, befintligt index",
       caption = "Skuggade fält indikerar en standardavvikelse") +
  facet_wrap(~SkolSDO) +
  theme_minimal(base_size = 11) +
  theme_rise() +
  theme(text = element_text(family = "Lato"))

Code
# psfTidscore <- df.index %>%
#   filter(Kön %in% c("Pojke", "Flicka")) %>%
#   filter(!is.na(SkolSDO)) %>%
#   group_by(ar, SkolSDO, Kön) %>%
#   summarise(PSFmean = round(mean(PSMscore100, na.rm = T), 0),
#             stdev = round(sd(PSMscore100, na.rm = T))) %>%
#   ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
#   geom_line(linewidth = 1) +
#   geom_point_interactive(aes(tooltip = PSFmean),
#                          size = 1.5) +
#   geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
#               alpha = 0.15, linetype = 0) +
#   geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
#              linetype = 2, color = "black", alpha = 0.8) +
#   geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
#              linetype = 3, color = "black", alpha = 0.8) +
#   scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
#   scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
#   scale_y_continuous("Medelvärde", limits = c(0,100)) +
#   labs(title = "Psykosomatiska besvär indexvärde",
#        caption = "Skuggade fält indikerar en standardavvikelse") +
#   facet_wrap(~SkolSDO) +
#   theme_minimal(base_size = 11) +
#   theme_rise()
# 
# ggiraph(ggobj = psfTidscore)

df.index %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(ar, SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PSMscore100, na.rm = T), 0),
            stdev = round(sd(PSMscore100, na.rm = T))) %>%
  ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
              alpha = 0.15, linetype = 0) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
             linetype = 2, color = "black", alpha = 0.8) +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
             linetype = 3, color = "black", alpha = 0.8) +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
  scale_y_continuous("Medelvärde", limits = c(0,100)) +
  labs(title = "Psykosomatiska besvär indexvärde",
              subtitle = "5 items, nytt index",
       caption = "Skuggade fält indikerar en standardavvikelse") +
  facet_wrap(~SkolSDO) +
  theme_minimal(base_size = 11) +
  theme_rise() +
  theme(text = element_text(family = "Lato"))

Code
df.index8 %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(ar, SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PS8score100, na.rm = T), 0),
            stdev = round(sd(PS8score100, na.rm = T))) %>%
  ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
              alpha = 0.15, linetype = 0) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
             linetype = 2, color = "black", alpha = 0.8) +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
             linetype = 3, color = "black", alpha = 0.8) +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
  scale_y_continuous("Medelvärde", limits = c(0,100)) +
  labs(title = "Psykiska/psykosomatiska besvär indexvärde",
              subtitle = "8 items, nytt index",
       caption = "Skuggade fält indikerar en standardavvikelse") +
  facet_wrap(~SkolSDO) +
  theme_minimal(base_size = 11) +
  theme_rise() +
  theme(text = element_text(family = "Lato"))

4.22.4 Fördelningar

Code
df.indexR %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PSFindex) %>% 
  summary()
    PSFindex     
 Min.   :  0.00  
 1st Qu.: 40.00  
 Median : 55.00  
 Mean   : 55.49  
 3rd Qu.: 75.00  
 Max.   :100.00  
Code
df.indexR %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PSFindex) %>% 
  pull() %>% 
  hist(., col = RISEprimGreen, main = "Psykisk ohälsa indexvärden")

Code
summary(df.indexR$PSFindex)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   40.00   60.00   57.54   75.00  100.00 
Code
df.index %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PSMscore100) %>% 
  summary()
  PSMscore100    
 Min.   :  0.00  
 1st Qu.: 41.00  
 Median : 51.00  
 Mean   : 50.29  
 3rd Qu.: 62.00  
 Max.   :100.00  
Code
df.index %>% 
  filter(ar == 2020) %>%  
  filter(ARSKURS == "Åk 9") %>% 
  select(PSMscore100) %>% 
  pull() %>% 
  hist(., col = RISEprimGreen, main = "Psykosomatiska besvär 5 item")

Code
df.index8 %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PS8score100) %>% 
  summary()
  PS8score100    
 Min.   :  0.00  
 1st Qu.: 43.00  
 Median : 50.00  
 Mean   : 50.87  
 3rd Qu.: 58.00  
 Max.   :100.00  
Code
df.index8 %>% 
  filter(ar == 2020) %>%
  filter(ARSKURS == "Åk 9") %>% 
  select(PS8score100) %>% 
  pull() %>% 
  hist(., col = RISEprimGreen, main = "Psykosomatiska besvär 8 item")

4.22.5 Violin- & boxplot jämförelse

Code
library(PrettyCols)
df.indexR %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutater = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSFindex, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "TealGreens") +
  labs(
    title = "Psykisk hälsa",
    subtitle = "Befintligt index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  )

Code
df.index %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutater = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSMscore100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "TealGreens") +
  labs(
    title = "Psykosomatiska besvär",
    subtitle = "Nytt 5-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  )

Code
df.index8 %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutater = factor(ar)) %>% 
  ggplot(aes(x = År, y = PS8score100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "TealGreens") +
  labs(
    title = "Psykiska/psykosomatiska besvär",
    subtitle = "Nytt 8-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  )

4.22.5.1 Fördelat på kön

Code
library(PrettyCols)
df.indexR %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutater = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSFindex, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "Rainbow") +
  labs(
    title = "Psykisk hälsa",
    subtitle = "Befintligt index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  ) + 
  facet_wrap(~Kön)

Code
df.index %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutater = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSMscore100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "Rainbow") +
  labs(
    title = "Psykosomatiska besvär",
    subtitle = "Nytt 5-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  ) +
  facet_wrap(~Kön)

Code
df.index8 %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutater = factor(ar)) %>% 
  ggplot(aes(x = År, y = PS8score100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "Rainbow") +
  labs(
    title = "Psykiska/psykosomatiska besvär",
    subtitle = "Nytt 8-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  ) + 
  facet_wrap(~Kön)

4.22.6 Tak-/golveffekter

Code
df.indexR %>% 
  group_by(PSFindex) %>% 
  reframe(n = n()) %>% 
  mutate(Procent = round(100*n/sum(n),1)) %>% 
  formattable(., table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 50%"')
PSFindex n Procent
0 577 0.7
5 865 1.0
10 1213 1.4
15 1654 1.9
20 2332 2.6
25 2937 3.3
30 3580 4.0
35 4307 4.9
40 5332 6.0
45 5841 6.6
50 6639 7.5
55 7077 8.0
60 7320 8.3
65 7190 8.1
70 7028 7.9
75 6367 7.2
80 5752 6.5
85 4339 4.9
90 3362 3.8
95 2271 2.6
100 2537 2.9
Code
df.index %>% 
  group_by(PSMscore100) %>% 
  reframe(n = n()) %>% 
  mutate(Procent = round(100*n/sum(n),1)) %>% 
  formattable(., table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 50%"')
PSMscore100 n Procent
0 1958 2.0
18 2503 2.5
26 3512 3.6
32 4807 4.9
37 5733 5.8
41 6804 6.9
44 8239 8.4
48 9286 9.4
51 9385 9.5
54 9501 9.6
58 9154 9.3
62 8324 8.5
66 7015 7.1
72 5374 5.5
81 3564 3.6
100 3307 3.4
Code
df.index8 %>% 
  group_by(PS8score100) %>% 
  reframe(n = n()) %>% 
  mutate(Procent = round(100*n/sum(n),1)) %>% 
  formattable(., table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 50%"')
PS8score100 n Procent
0 338 0.4
15 572 0.6
22 876 0.9
27 1275 1.3
31 1753 1.8
34 2293 2.4
36 2832 3.0
39 3543 3.7
41 3954 4.1
43 4712 4.9
45 5401 5.6
47 5822 6.1
48 6385 6.7
50 6624 6.9
52 6834 7.1
54 6938 7.2
56 6521 6.8
58 6218 6.5
61 5801 6.0
64 5021 5.2
67 4096 4.3
71 3373 3.5
76 2297 2.4
84 1467 1.5
100 1008 1.1

4.23 Programvara som använts för analyserna

Code
pkgs <- cite_packages(cite.tidyverse = TRUE, 
                      output = "table",
                      bib.file = "grateful-refs.bib",
                      include.RStudio = TRUE)
formattable(pkgs, 
            table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 80%"')
Package Version Citation
arrow 10.0.0 Richardson et al. (2022)
base 4.2.2 R Core Team (2022)
car 3.1.1 Fox and Weisberg (2019)
catR 3.17 Magis and Raîche (2012); Magis and Barrada (2017)
colorspace 2.1.0 Zeileis, Hornik, and Murrell (2009); Stauffer et al. (2009); Zeileis et al. (2020)
cowplot 1.1.1 Wilke (2020)
eRm 1.0.2 Mair and Hatzinger (2007b); Mair and Hatzinger (2007a); Hatzinger and Rusch (2009); Rusch, Maier, and Hatzinger (2013); Koller, Maier, and Hatzinger (2015); Debelak and Koller (2019); Mair, Hatzinger, and Maier (2021)
extrafont 0.18 Chang (2022)
foreach 1.5.2 Microsoft and Weston (2022)
formattable 0.2.1 Ren and Russell (2021)
furrr 0.3.1 Vaughan and Dancho (2022)
ggchicklet 0.6.0 Rudis (2022)
ggdist 3.2.0 Kay (2022)
gghalves 0.1.4 Tiedemann (2022)
ggpp 0.4.5 Aphalo (2022)
ggrepel 0.9.3 Slowikowski (2023)
glue 1.6.2 Hester and Bryan (2022)
grateful 0.1.11 Rodríguez-Sánchez, Jackson, and Hutchins (2022)
kableExtra 1.3.4 Zhu (2021)
knitr 1.42 Xie (2014); Xie (2015); Xie (2023)
matrixStats 0.63.0 Bengtsson (2022)
mirt 1.37.1 Chalmers (2012)
PrettyCols 1.0.1 Rennie (2023)
psych 2.2.9 Revelle (2022)
psychotree 0.16.0 Trepte and Verbeet (2010); Strobl, Wickelmaier, and Zeileis (2011); Strobl, Kopf, and Zeileis (2015); Komboz, Zeileis, and Strobl (2018); Wickelmaier and Zeileis (2018)
reshape 0.8.9 Wickham (2007)
RISEkbmRasch 0.1.11.1 Johansson (2023)
rmarkdown 2.20 Xie, Allaire, and Grolemund (2018); Xie, Dervieux, and Riederer (2020); Allaire et al. (2023)
scales 1.2.1 Wickham and Seidel (2022)
tidyverse 2.0.0 Wickham et al. (2019)

4.24 Referenser

Allaire, JJ, Yihui Xie, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang, and Richard Iannone. 2023. Rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.
Aphalo, Pedro J. 2022. Ggpp: Grammar Extensions to ’Ggplot2’. https://CRAN.R-project.org/package=ggpp.
Bengtsson, Henrik. 2022. matrixStats: Functions That Apply to Rows and Columns of Matrices (and to Vectors). https://CRAN.R-project.org/package=matrixStats.
Chalmers, R. Philip. 2012. mirt: A Multidimensional Item Response Theory Package for the R Environment.” Journal of Statistical Software 48 (6): 1–29. https://doi.org/10.18637/jss.v048.i06.
Chang, Winston. 2022. Extrafont: Tools for Using Fonts. https://CRAN.R-project.org/package=extrafont.
Debelak, Rudolf, and Ingrid Koller. 2019. Testing the Local Independence Assumption of the Rasch Model With Q3-Based Nonparametric Model Tests.” Applied Psychological Measurement. https://doi.org/10.1177/0146621619835501.
Fox, John, and Sanford Weisberg. 2019. An R Companion to Applied Regression. Third. Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.
Hatzinger, Reinhold, and Thomas Rusch. 2009. IRT models with relaxed assumptions in eRm: A manual-like instruction.” Psychology Science Quarterly 51.
Hester, Jim, and Jennifer Bryan. 2022. Glue: Interpreted String Literals. https://CRAN.R-project.org/package=glue.
Johansson, Magnus. 2023. RISEkbmRasch: Psychometric Analysis in r with Rasch Measurement Theory. https://github.com/pgmj/RISEkbmRasch.
Kay, Matthew. 2022. ggdist: Visualizations of Distributions and Uncertainty. https://doi.org/10.5281/zenodo.3879620.
Koller, Ingrid, Marco Johannes Maier, and Reinhold Hatzinger. 2015. An Empirical Power Analysis of Quasi-Exact Tests for the Rasch Model: Measurement Invariance in Small Samples.” Methodology 11. https://doi.org/10.1027/1614-2241/a000090.
Komboz, Basil, Achim Zeileis, and Carolin Strobl. 2018. “Tree-Based Global Model Tests for Polytomous Rasch Models.” Educational and Psychological Measurement 78 (1): 128–66. https://doi.org/10.1177/0013164416664394.
Magis, David, and Juan Ramon Barrada. 2017. “Computerized Adaptive Testing with R: Recent Updates of the Package catR.” Journal of Statistical Software, Code Snippets 76 (1): 1–19. https://doi.org/10.18637/jss.v076.c01.
Magis, David, and Gilles Raîche. 2012. “Random Generation of Response Patterns Under Computerized Adaptive Testing with the R Package catR.” Journal of Statistical Software 48 (8): 1–31. https://doi.org/10.18637/jss.v048.i08.
Mair, Patrick, and Reinhold Hatzinger. 2007a. CML based estimation of extended Rasch models with the eRm package in R.” Psychology Science 49.
———. 2007b. Extended Rasch modeling: The eRm package for the application of IRT models in R.” Journal of Statistical Software 20. https://www.jstatsoft.org/v20/i09.
Mair, Patrick, Reinhold Hatzinger, and Marco Johannes Maier. 2021. eRm: Extended Rasch Modeling. https://cran.r-project.org/package=eRm.
Microsoft, and Steve Weston. 2022. Foreach: Provides Foreach Looping Construct. https://CRAN.R-project.org/package=foreach.
R Core Team. 2022. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Ren, Kun, and Kenton Russell. 2021. Formattable: Create ’Formattable’ Data Structures. https://CRAN.R-project.org/package=formattable.
Rennie, Nicola. 2023. PrettyCols: Pretty Colour Palettes. https://CRAN.R-project.org/package=PrettyCols.
Revelle, William. 2022. Psych: Procedures for Psychological, Psychometric, and Personality Research. Evanston, Illinois: Northwestern University. https://CRAN.R-project.org/package=psych.
Richardson, Neal, Ian Cook, Nic Crane, Dewey Dunnington, Romain François, Jonathan Keane, Dragoș Moldovan-Grünfeld, Jeroen Ooms, and Apache Arrow. 2022. Arrow: Integration to ’Apache’ ’Arrow’. https://CRAN.R-project.org/package=arrow.
Rodríguez-Sánchez, Francisco, Connor P. Jackson, and Shaurita D. Hutchins. 2022. Grateful: Facilitate Citation of r Packages. https://github.com/Pakillo/grateful.
Rudis, Bob. 2022. Ggchicklet: Create ’Chicklet’ (Rounded Segmented Column) Charts. https://git.rud.is/hrbrmstr/ggchicklet.
Rusch, Thomas, Marco Johannes Maier, and Reinhold Hatzinger. 2013. Linear logistic models with relaxed assumptions in R.” In Algorithms from and for Nature and Life, edited by Berthold Lausen, Dirk van den Poel, and Alfred Ultsch. Studies in Classification, Data Analysis, and Knowledge Organization. New York: Springer. https://doi.org/10.1007/978-3-319-00035-0_34.
Slowikowski, Kamil. 2023. Ggrepel: Automatically Position Non-Overlapping Text Labels with ’Ggplot2’. https://CRAN.R-project.org/package=ggrepel.
Stauffer, Reto, Georg J. Mayr, Markus Dabernig, and Achim Zeileis. 2009. “Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations.” Bulletin of the American Meteorological Society 96 (2): 203–16. https://doi.org/10.1175/BAMS-D-13-00155.1.
Strobl, Carolin, Julia Kopf, and Achim Zeileis. 2015. “Rasch Trees: A New Method for Detecting Differential Item Functioning in the Rasch Model.” Psychometrika 80 (2): 289–316. https://doi.org/10.1007/s11336-013-9388-3.
Strobl, Carolin, Florian Wickelmaier, and Achim Zeileis. 2011. “Accounting for Individual Differences in Bradley-Terry Models by Means of Recursive Partitioning.” Journal of Educational and Behavioral Statistics 36 (2): 135–53. https://doi.org/10.3102/1076998609359791.
Tiedemann, Frederik. 2022. Gghalves: Compose Half-Half Plots Using Your Favourite Geoms. https://CRAN.R-project.org/package=gghalves.
Trepte, Sabine, and Markus Verbeet, eds. 2010. Allgemeinbildung in Deutschland – Erkenntnisse Aus Dem SPIEGEL Studentenpisa-Test. Wiesbaden: VS Verlag.
Vaughan, Davis, and Matt Dancho. 2022. Furrr: Apply Mapping Functions in Parallel Using Futures. https://CRAN.R-project.org/package=furrr.
Wickelmaier, Florian, and Achim Zeileis. 2018. “Using Recursive Partitioning to Account for Parameter Heterogeneity in Multinomial Processing Tree Models.” Behavior Research Methods 50 (3): 1217–33. https://doi.org/10.3758/s13428-017-0937-z.
Wickham, Hadley. 2007. “Reshaping Data with the Reshape Package.” Journal of Statistical Software 21 (12). https://www.jstatsoft.org/v21/i12/.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Wickham, Hadley, and Dana Seidel. 2022. Scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.
Wilke, Claus O. 2020. Cowplot: Streamlined Plot Theme and Plot Annotations for ’Ggplot2’. https://CRAN.R-project.org/package=cowplot.
Xie, Yihui. 2014. “Knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.
———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.
———. 2023. Knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.
Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.
Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.
Zeileis, Achim, Jason C. Fisher, Kurt Hornik, Ross Ihaka, Claire D. McWhite, Paul Murrell, Reto Stauffer, and Claus O. Wilke. 2020. colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes.” Journal of Statistical Software 96 (1): 1–49. https://doi.org/10.18637/jss.v096.i01.
Zeileis, Achim, Kurt Hornik, and Paul Murrell. 2009. “Escaping RGBland: Selecting Colors for Statistical Graphics.” Computational Statistics & Data Analysis 53 (9): 3259–70. https://doi.org/10.1016/j.csda.2008.11.033.
Zhu, Hao. 2021. kableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.